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Numerical schemes used for the integration of complex flow simulations should provide accurate solutions for the long time in- 
tegrations these flows require. To this end, the performance of various high-order accurate numerical schemes is investigated for 
direct numerical simulations (DNS) of homogeneous isotropic two-dimensional decaying turbulent flows. The numerical accuracy 
of compact difference, explicit central difference, Arakawa, and dispersion-relation-preserving schemes are analyzed and compared 
with the Fourier- Galerkin pseudospectral scheme. In addition, several explicit Runge-Kutta schemes for time integration are inves- 
tigated. We demonstrate that the centered schemes suffer from spurious Nyquist signals that are generated almost instantaneously 
and propagate into much of the field when the numerical resolution is insufficient. We further show that the order of the scheme be- 
comes increasingly important for increasing cell Reynolds number. Surprisingly, the sixth-order schemes are found to be in perfect 
agreement with the pseudospectral method. Considerable reduction in computational time compared to the pseudospectral method 
is also reported in favor of the finite difference schemes. Among the fourth-order schemes, the compact scheme provides better 
accuracy than the others for fully resolved computations. The fourth-order Arakawa scheme provides more accurate results for 
under-resolved computations, however, due to its conservation properties. Our results show that, contrary to conventional wisdom, 
difference methods demonstrate superior performance in terms of accuracy and efficiency for fully resolved DNS computations of 
the complex flows considered here. For under-resolved simulations, however, the choice of difference method should be made with 
care. 

Keywords: Incompressible flows, high-order accurate schemes, finite difference methods, Fourier- Galerkin pseudospectral 
method, double shear layer problem, two-dimensional decaying turbulence 
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1. Introduction 

The physics of two-dimensional turbulence have been elu- 
cidated substantially during the past decades by theoretical 
models, intensive numerical investigations, and dedicated soap 
film experiments [ 1 ]. Two-dimensional turbulence research ef- 
forts have applicability in geophysics, astronomy and plasma 
physics, in which numerical experiments play a large role. One 
of the most important reasons for studying two-dimensional tur- 
bulence is to improve our understanding of geophysical flows 
in the atmosphere and ocean EHEJ. We may also find two- 
dimensional flows in a wide variety of situations such as flows 
in rapidly rotating systems and flows in a fluid film on top of 
the surface of another fluid or a rigid object |9). 

Two-dimensional turbulence behaves in a profoundly differ- 
ent way from three-dimensional turbulence due to different en- 
ergy cascade behavior, and follows the Kraichnan-Batchelor- 
Leith (KBL) theory fT0lfT2l . In three-dimensional turbulence, 
energy is transferred forward, from large scales to smaller 
scales, via vortex stretching. In two dimensions that mechanism 
is absent, and it turns out that under most forcing and dissipa- 
tion conditions energy will be transferred from smaller scales 



* Corresponding author. 

Email address: omersan@vt . edu (Omer San) 



to larger scales. This is largely because of another quadratic 
invariant, the potential enstrophy, defined as the integral of the 
square of the potential vorticity. Despite the apparent simplicity 
in dealing with two rather than three spatial dimensions, two- 
dimensional turbulence is possibly richer in its dynamics than 
three-dimensional turbulence due to its conservation properties, 
such as its inverse energy and forward enstrophy cascading 
mechanisms. Danilov and Gurarie fT3l and Tabeling lfl4l re- 
viewed both theoretical and experimental two-dimensional tur- 
bulence studies along with extensions into geophysical flow set- 
tings. More recent reviews on two-dimensional turbulence are 
also provided by Clercx and van Heijst 031 and Boffetta and 
Ecke 1 16 ]. Recent studies in two-dimensional turbulence, both 
forced (stationary) turbulence fT7H20l and unforced (decaying) 
turbulence l2TH23ll provide high resolution computational con- 
firmation of the KBL theory. 

Simulation of turbulent and other convection-dominated un- 
steady flows using direct numerical simulation (DNS) requires 
a numerical method that properly resolves all the multiscale 
flow structures |24|. Since high accuracy is crucial in numer- 
ical simulation of complex flows with multiscale structures, 
such as the unsteady evolution of a turbulent flow field, most 
two-dimensional turbulence studies have been performed us- 
ing pseudospectral methods based on fast Fourier transform 
(FFT) algorithms [16,21]. Simulations performed by the lattice 
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Boltzmann method (LBM) have been also presented for two- 
dimensional decaying turbulence |25 ]. Pseudospectral methods 
are highly accurate but mostly limited to ideal geometries such 
as rectangular or circular domains. Discretization methods such 
as finite difference, finite element, or finite volume methods are 
often preferred in more realistic problems. Finite difference 
methods offer an attractive alternative to spectral methods in 
the direct and large eddy simulations (LES) of turbulence pro- 
viding reasonable accuracy coupled with relative ease of im- 
plementation in simple and complex flow geometries 126^291 . 
Computational algorithms developed in the past were mainly 
designed for solving large-scale fluid dynamics problems us- 
ing second-order spatial accuracy [30-32]. These algorithms 
usually have rather significant dispersion errors and if they are 
not centered schemes they also have large dissipation errors, 
making it hard to accurately compute fine structures in the flow 
field using them (33). There are two ways to improve the reso- 
lution of these methods; one is to refine the grid and the other is 
to construct a high-order accurate scheme. Our approach here 
is to test and evaluate different high-order formulations for in- 
stantaneous and statistical properties of two-dimensional tur- 
bulence and compare their accuracy and efficiency with those 
of the pseudospectral and the second-order schemes. Further- 
more, it has been shown by Kravchenko and Moin l34l that the 
subgrid- scale models in LES are effective only if central dis- 
cretization of order higher than two is employed. With this in 
mind, we will investigate the behavior of four different families 
of high-order accurate finite difference methods in the decay of 
two-dimensional isotropic turbulence. 

High-order finite difference schemes can be formulated to 
reduce the truncation errors associated with the difference ap- 
proximations. A straightforward Taylor series expansion of a 
pointwise discretization under certain assumptions results in 
a family of the explicit difference (ED) schemes. The com- 
pact difference (CD) schemes feature high-order accuracy with 
smaller stencils and smaller truncation errors than the ED 
schemes, and have been employed as an alternative to spec- 
tral methods in simulations of turbulence with great flexibility 
l35l . On the other hand, increasing the stencil size allows us 
to optimize the weight coefficients in the difference equation. 
This strategy leads to the dispersion-relation-preserving (DRP) 
schemes l36lL which have been used mostly in acoustics. An- 
other strategy to construct a numerical scheme is based on the 
conservation properties of the discrete form of the equations. 
Arakawa l37l suggested that the conservation of energy, enstro- 
phy, and skew- symmetry is sufficient to avoid computational in- 
stabilities stemming from nonlinear interactions. The conserva- 
tion and stability properties of the Arakawa scheme were inves- 
tigated by Lilly 0T1 by means of spectral analysis along with 
several first and second-order time integration methods. In the 
present work, we test several Runge-Kutta methods for time in- 
tegration, although the primary goal here is to analyze the accu- 
racy of these high-order accurate spatial differencing methods 
for the long-term evolution of complex two-dimensional turbu- 
lent flows. For finite difference schemes, the combination of 
differentiation errors and nonlinear truncation and aliasing er- 
rors, which usually manifest themselves in the high wavenum- 



bers of the resolved scales, determines the overall error at the 
small scales. Looking at the accuracy of the whole solution pro- 
cedure we also investigate the resolution requirements for these 
finite difference families, the effects of the order of the schemes, 
and the importance of the global conservation properties. 

The paper is organized as follows: the mathematical for- 
mulation of the problem is given in Section [2] The numer- 
ical methods are presented in Section [3] with descriptions of 
high-order accurate spatial discretization schemes, temporal 
discretization algorithms, and an efficient fast Poisson solver 
algorithm. These schemes are validated in Section [4] by sim- 
ulating the Taylor-Green decaying vortex benchmark problem 
for the unsteady incompressible Navier-Stokes equations. The 
effective accuracies of these methods are also provided in this 
section, and are confirmed to be the theoretical accuracies of 
the schemes. Section [5] presents a careful numerical investiga- 
tion of their performance for a challenging benchmark prob- 
lem which consists of strong shear layers. The results for two- 
dimensional isotropic homogeneous decaying turbulence are 
provided in Section [6] The behavior of these nine different spa- 
tial schemes are tested in terms of accuracy and efficiency. The 
effects of several explicit Runge-Kutta time advancement tech- 
niques on the whole solution procedure are also analyzed. In 
addition, the Reynolds number (Re) dependency of the turbu- 
lence statistics is illustrated in this section. Final conclusions 
and some comments on the performance of these schemes are 
drawn in section 

2. Mathematical model 

The governing equations for two-dimensional incompress- 
ible flows can be written in a dimensionless form of the 
vorticity- stream function formulation as 

doj ^ dif/ doj dif/ dco 1 ^d 2 oj ^ d 2 oj^ 
dt dy dx dx dy Re dx 2 dy 2 

along with the kinematic relationship between vorticity and 
stream function according to a Poisson equation, which is given 
as 

From a computational point of view, this formulation has sev- 
eral advantages over the primitive variable formulation. It elim- 
inates pressure from the Navier-Stokes equations and hence 
has no corresponding odd-even decoupling between the pres- 
sure and velocity components, as well as projection inaccura- 
cies usually observed in fractional step approaches (38). There- 
fore, the usage of a collocated grid does not produce any 
spurious modes in the vorticity- stream function formulation. 
The vorticity- stream function formulation automatically satis- 
fies the divergence-free condition and allows one to reduce the 
number of equations to be solved. 

The main objective of our work is to test and evalu- 
ate different frameworks for high-order accurate finite differ- 
ence schemes and compare them with a spectrally accurate 
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pseudospectral method for two-dimensional isotropic turbulent 
flows. In fact, to be able to compare the numerical schemes 
more precisely we restricted ourselves to periodic boundary 
conditions and a uniform Cartesian grid. Consequently, we 
eliminated errors coming from the mesh non-uniformities and 
inconsistent boundary schemes. It should also be noted that 
using the vorticity-stream function formulation on a collocated 
grid provides us with an ideal computational setting in which to 
test the characteristics of the numerical schemes by eliminating 
any possible errors coming from projection inaccuracies. 

3. Numerical methods 

The objective of the present work is to test and evaluate 
different frameworks for high-order accurate finite difference 
schemes and compare them with a spectrally accurate pseu- 
dospectral method for two-dimensional isotropic turbulence 
flows. We will also compare them with their classical second- 
order accurate formulations in order to analyze the effects of 
order of accuracy. In the text that follows, the high-order finite 
difference approximations for spatial discretization, the pseu- 
dospectral method, an efficient procedure for the Poisson equa- 
tion, and the time integration methods that we utilize for time 
advancement are given briefly for completeness. 



which gives rise to an a-family of tridiagonal schemes with 
a = | (or + 2), and b = \{Aa - 1). Here, a = leads to the 
explicit non-compact fourth-order scheme for the first deriva- 
tive. A classical compact fourth-order scheme, which is also 
known as the Pade scheme, is obtained by setting a - 1/4. The 
truncation error in the Eq. (|9j) is ^ (3a - l)/z 4 / (5) . Therefore, a 
sixth-order compact scheme is obtained by choosing a = 1/3. 
The second derivative compact centered scheme is given by 



af^+ff' + afll^a 



f M - 2fi + fi-i 



h 2 

f i+2 ~ 2f, + fi-2 



Ah 2 



(10) 



|(1 - a), and b = ±(10a - 1). For a = 1/10 the 



where a 

classical fourth-order Pade scheme is obtained. The truncation 
error in the Eq. (flOb is -^(Ua - 2)h A f {6 \ Therefore, a sixth- 
order compact scneme is also obtained by choosing a = 2/1 1. 

These compact schemes for the first and second derivatives 
show better spectral accuracy then their explicit counterparts. 
However, they involve solving tridiagonal matrix equations and 
are therefore effectively non-local. This requires special care 
for designing parallel computations. However, in this study, we 
focus on a single domain without domain decomposition, and 
so the compact schemes can be adopted easily. 



3. 1 . Explicit difference ( ED ) scheme 

For any scalar pointwise value of / the classical centered fi- 
nite difference schemes for the first derivative up the the sixth 
order accuracy are given by |39l 



fl =^(-fi-i+f i+ i) + 0(h 2 ) 

f'i = ^(fi-2 ~ 8/,-i + Sf i+1 - f i+2 ) + 0(h 4 ) 

f!= g^(-^-3+9^_ 2 -45/-_i 

+45f M -9f i+2 +f M ) + 0(h 6 ) 
Similarly, the second derivative approximations become 

f l " = ^(fi-i-2f i + f i+l ) + 0(h 2 ) 

f'i = J^(~f<-2 + 16 ^-1 - 3 °f' 



+I6f i+l - f i+2 ) + 0(h 4 ) 
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f'i = ~ 2 ' 7 fi- 2 + 270/i-i - 490/i 

+210f M - 27f i+2 + 2f i+3 ) + 0(h 6 ) 
where h is the step size in the derivative direction. 



(3) 



(4) 



(5) 



(6) 



(7) 



(8) 



3.2. Compact difference ( CD) scheme 

Apart from the explicit centered schemes, with the cost of 
a tridiagonal system solution (i.e., by using the well-known 
Thomas algorithm), the first derivative can be computed accord- 
ingly 



<xfi_i +&+ af M = a + f 



3.3. Dispersion-relation-preserving (DRP) scheme 

An optimized dispersion-preserving scheme was proposed in 
l36l . The scheme increases the stencil size, but improves the 
resolution in wavespace. The finite difference approximation 
for the first derivative is 



1 M 
K = hTj a ' fi+ > 



(11) 



j=-N 



Fourth-order accuracy is obtained for N = M = 3 (centered 
DRP scheme) and the optimized coefficients are a$ = 0, a\ = 
-a- x = 0.79926643, a 2 = -a. 2 = -0.18941314, and a 3 = 
-a-3 = 0.02651995. The viscous terms can be computed by 
using the explicit difference schemes. 

3.4. Arakawa scheme 

Arakawa 1 37 ] suggested that the conservation of energy, en- 
strophy and skew symmetry is sufficient to avoid computational 
instabilities arising from nonlinear interactions. The nonlinear 
terms in Eq. ([T]) were defined as the Jacobian 

dif/ doj dif/ dco 
dy dx dx dy 

The second-order Arakawa scheme for the Jacobian is 



J(0), iff) 



(12) 



2h 



Ah 



(9) 



JK^iff) = ^(Ji(a),iff) + JiM) + U^iff)) (13) 
where the discrete parts of the Jacobian 

Ji(u,iff) = y^[(a> i+ ij - a)i-ij)(iffij +1 - iffij-i) 

- (a)ij + i - a)ij-i)(iff i+ ij - iffi-ij)] (14) 



J 2 (oj,iff) = -——[oj i+ ij(i// i+ ij + i -ifs i+ ij-i) 



4h x h y 



- OJi-ijilffi-ij+i - l/fi-lj-l) 



(15) 



J 3 (0), iff) 



1 



Ah x hy 



[a) i+1J+ i(ifs iJ+ i -iffi+ij) 



- ^-ij-iC^-ij - iffij-i) 

- (Oi-ij+iifcj+i -iffi-lj) 

+ a)i + ij-i(!ff i+ ij - iffij-i)] 



(16) 



The fourth-order accurate Arakawa discretization of the Jaco- 
bian becomes 

J n (a), iff) = ^(Ma>, iff) + J 5 (o), + Je(o), <A)) (17) 



where 

J 4 (co, iff) 



8h x h y 



[(0Ji + ij + i - (i)i-ij-i)(lffi-ij+i - iffi+ij-i) 



- ((Oi-ij+i - o)i+ij-i)(iffi + ij+i - iffi-ij-i)] (18) 



J 5 (0J,lff) = ——-[a) i+ ij + i(lffij + 2 - lffi+2,j) 

on x n y 

- QJi-ij-i(lffi- 2 J - iffij-l) 

- (i)i-ij + i(lff Uj+ 2 ~ lffi-2,j) 
+ 0) i+ i j-i(lff i+2J - lffi,j- 2 )] 



(19) 



J 6 (oj, if/) 



1 



8h x h y 



[CL)i + 2,j(lffi+l,j+\ -lffi + l,j-l) 



- COi-2,j(lffi-l,j+l ~ fc-lj-l) 
~ MiJ+lifti+lJ+l ~ fti-lj+l) 

+ a)ij-2(iffi+ij-i - if/i-ij-i)] 



(20) 



Arakawa showed that /// conserves enstrophy and energy and 
the following Jacobian 



J(oj, iff) = Ufa, iff) - J n (a), iff) + OQT) 



(21) 



has fourth-order accuracy. The viscous terms can be discretized 
with the 4th-order explicit difference scheme given in the previ- 
ous section. This scheme was used to compute two-dimensional 
isotropic turbulence by Orlandi [ 24 1. 

3.5. Fourier- Gale rkin pseudo spectral method 

Fourier series expansion based methods are often used for 
solving problems with periodic boundary conditions. One of 
the most accurate methods for solving the Navier-Stokes equa- 
tions in periodic domains is the pseudospectral method, which 
exploits fast Fourier transform (FFT) algorithms, resulting in 
spectral accuracy l40l . It should be noted that the discrete 
Fourier transform is 0(N^Ny), which becomes very expensive 
for larger resolutions. On the other hand, the discrete Fourier 
transform can, in fact, be computed in 0(N x N y log(N x N y )) oper- 
ations with well-known FFT algorithms BTll . Given a set of N x 



and N y uniformly distributed points on the interval [0, L x ]x [0, 
Ly], the forward Fourier transformation of any discrete function 



u is 



Nx_i Ny _ 



2_j u m, n e Lx 



yj) 



(22) 



and its inverse transform to find the Fourier coefficients is 



11 



N x -\N y -\ 

wn -—^y Y Uije-^^ (23) 



where i = V-T and i, and j represent indices for physical space, 
and m, and n are the Fourier space indices. The discrete grid 
coordinates for i = 0, 1, ..N x and j = 0, 1, ..N y are given by 



N v ' 



(24) 



Since the domain is periodic, xo = x^ x and yo = y^ y . Defining 
the wave numbers 



2nm 2nn 



(25) 



Now, we can easily perform differentiation in both directions. 
The first and second order differentiation of any function u in 
discrete domain becomes 



ditij 
dx 



+k y yj) 



c ru-, 



1 o 



dx 2 



r = Z Z ~ U ^-£ X Y 



JL{k x Xi+kyyj) 



(26) 



(27) 



and similarly for the y direction derivatives. By transforming 
Eq. ([]]> and Eq. (|2]i to Fourier space, the governing equations 
become 



^T+N= ^[(-k 2 x -kl)cd m , n ] 
at Re y 



(28) 
(29) 



where N represents the nonlinear terms (Fourier transform of 
the Jacobian) and needs to be computed according to the fol- 
lowing convolution 

N = (ik y iff m ^) o (ik x a) m , n ) - (ik x iff m , n ) o (ik y (b m , n ) (30) 

The convolution sum for these nonlinear terms is actually com- 
puted in the physical domain, with the help of the convolu- 
tion theorem. One of the basic techniques for removing the 
aliasing error for the nonlinear terms is the 3/2-rule which is 
known as the padding or truncation algorithm. The key to this 
dealiasing algorithm is the use of the discrete transform of M x 
and My rather than N x and N y points, so that M x > 3N x /2 
and My > 3N y /2. For example, let's denote a m ^ n and b m ^ n , 



the Fourier coefficients of two multiplicand terms, as being in 
N, and we would like to compute the Fourier coefficients of 
the multiplication of these two terms, c m>n = a m ^ n o b m , n . In 
the pseudospectral (PS) method, we can compute this term as 
c m ,n - F~ l (F(a m , n )F(b m , n )), where F and F~ l represent the 
forward and inverse transforms, respectively. First we need to 
use the forward transformation using M x and M y wavenumbers 
such that 



KkxXi+kyYj) 



A - V V ~ 1 



(31) 



(32) 



where X t = ^ , Yj = ^ . The Fourier coefficients are padded 
with zeros for the additional wavenumbers such as 



a m , n \m\ < N x , \n\ < N y 
otherwise 







\m\ < N x , \n\ < N y 
otherwise 



(33) 



(34) 



Multiplication is performed in physical space but with padded 
higher wavenumbers in the following way 



Cij - AijBij. 



(35) 



Then, the inverse transform gives the dealiased Fourier coeffi- 
cients as 



r 



1 



M X -\ My~] 



z z 



l(k x X i+ k Y Yj) 



(36) 



in which we are only interested in C m?w for \m\ < N x , \n\ < N y 
(i.e., c m , n - C mrn for \m\ < N x , \n\ < N y ). Although in practice 
dealiasing need only be performed with M x > 3N x /2, we have 
used M x = 2N X here since the FFT solver we are using requires 
resolutions that are powers of two. In order to make a fair com- 
parison of CPU times, however, we have chosen to present the 
dealisased computational CPU times an appropriately interpo- 
lated factor of the actual time between N x and M x . In this study, 
the semi-discrete vorticity transport equation in Fourier space, 



Eq. (28 ), is solved explicitly by third-order or forth-order accu- 
rate Runge-Kutta time marching schemes which will be given 
later. 

3.6. Poisson solver 

In order to obtain the streamfunction from the vorticity field, 
Eq. ^ needs to be solved. In general, Gauss-Seidel or suc- 
cessive over relaxation (SOR) types of iterative algorithms for 
solving the Poisson equation are of 0(N 2 ) where N is the to- 
tal number of grid points (i.e., TV = N x N y for two-dimensional 
problems). It is not feasible to use these types of iterative Pois- 
son solvers for high resolution (and therefore high Reynolds 
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Figure 1 : Efficiency of Poisson solvers for a two-dimensional square domain. 



number) computations along with long time integration. In 
order to accelerate these solvers, very successful multi-grid 
(MG) algorithms have been developed that reduce the compu- 
tational effort to 0(CmgN) where Cmg is a proportionality con- 
stant 1 42). On the other hand, for certain ideal problems on 
equally- spaced grids, fast Fourier transform (FFT)-based fast 
Poisson solvers (FPS) can be used to solve elliptic problems di- 
rectly with an 0(C FP sNlog(N)) operational cost, and they are 
presently the fastest available algorithms (CfpslogiN) < Cmg 
in relevant resolutions) for solving Poisson equations BT1I431 . 

A preliminary study is performed to test the computational 
efficiencies of different Poisson solvers for a square domain 
with equidistant grid spacing. The computational time (CPU 
time) for solving just one Poisson equation is illustrated in 
Fig. [T] This preliminary comparison shows that the FPS is the 
most efficient solver for Cartesian grid problems. We take ad- 
vantage of the simple rectangular shape of our domain and uti- 
lize the FPS for all the simulations in this study. The basic 
three- step procedure to find if/ from known oj values at any time 
(also for the inner steps of the Runge-Kutta time stepping algo- 
rithms) has the following form. 

1 . Use forward FFT for a> to find Fourier coefficients accord- 



ing to Eq. ( |22| ). 

2. Compute Fourier coefficients of stream function from 



Eq. (29) 



(37) 



3. Use inverse FFT for if/ values according to Eq. ( [23] ). 

3. 7. Time integration algorithms 

Steady calculations in fluid flows have been mostly handled 
using implicit time advancement algorithms and correspond- 
ing large time steps. It is tempting to adopt a similar practice 
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in computations of turbulent flows. Unfortunately, the require- 
ment of time accuracy over a wide range of scales does not 
permit very large time steps in such calculations [44 1. The 
use of large time steps implies that the small scales can have 
large errors, which can corrupt the solutions. Therefore, the 
families of fourth-order B31 and third-order l46l Runge-Kutta 
time integration schemes have been extensively used in turbu- 
lence simulations with their certain numerical stability require- 
ments [47 1 . In fluid dynamics, the numerical stability of such 
explicit algorithms is usually addressed by using the connec- 
tive Courant-Friedrichs-Lewy (CFL) or viscous Neumann sta- 
bility conditions |48l . Usually the fourth-order Runge-Kutta 
schemes have a greater accuracy for a given time step and a 
larger allowable stability domain than the third-order Runge- 
Kutta schemes. To be practical, an explicit algorithm for turbu- 
lence computations must provide long-time temporal accuracy 
while requiring modest temporary storage resulting in the de- 
velopment of a low-storage Runge-Kutta family for demanding 
computations. 

Semi-discrete ordinary differential equations (ODEs) are ob- 
tained after spatial discretization of the partial differential equa- 
tions. To implement the Runge-Kutta schemes, we cast the gov- 
erning vorticity transport equation in the following form 



Table 1: Coefficients for various low-storage third-order Runge-Kutta methods. 



scheme 


®2 


a 3 


Pi 


Pi 


P3 


symmetric 


2 
3 


-1 


1 

3 


1 


1 

2 


predictor/ corrector 


1 

4 


4 

3 


1 

2 


2 
3 


1 

2 


inhomogeneous 


17 

32 


32 
27 


1 

4 


8 
9 


3 
4 


Williamson (46J 


5 
9 


153 
128 


1 

3 


15 
16 


8 

15 



at time t n+l = f + At we set in Eq. ([41]) oj (0) = of and, af- 
ter the last step, of +l = 



.i(3) 



In order to be able to calcu- 



late the first step, i = 1, for which no g (U) exists, we require 
ol\ = 0. The corresponding coefficients for different types of 
constructed schemes are listed in Table[T] In practice all of them 
are approximately equally good as advocated by Brandenburg 
l5TTl . Although our main objective is to evaluate the perfor- 
mance of the various high-order spatial discretization schemes, 
we will perform a comparative study on these explicit Runge- 
Kutta schemes presented above to test their long-term charac- 
teristics for two-dimensional turbulence. 



doj 
dt 



: £(w, iff) 



(38) 4. Validation problem: Taylor-Green vortex 



where £(cj, if/) is the discrete operator of spatial derivatives in- 
cluding nonlinear convective terms and linear diffusive terms, 
and if/ is obtained from the Poisson equation. A fourth-order 
Runge-Kutta scheme can be written in the following form l49l 





= cf 




of" 


= co n 


At , , 

+ y£(o>V) 




= cf 


+ Af£(w 2 ,^ 2 ) 




1 

= 3 ( 





^£( W (3) ,.A (3) )(39) 
6 

The optimal third-order accurate TVD Runge-Kutta scheme is 
given as l50ll 



aP = co n + At£{co\if/ n ) 

-of + -co ( ' + - 
4 4 4 



(2) - Z^n + ^ (1) + ^£(^(1)^(1)) 



^ = * n + -J 2) + - At£(J 2) ,<A (2) ). 
3 3 3 r 



(40) 



Alternatively, memory-effective 2A/-storage schemes that re- 
quire only two sets of variables to be held in memory can be 
used. The iterative form of the 2yV-storage schemes can be writ- 
ten as 



g (0 = ^-l) +A ,£(^-l),^-l)) 



(41) 



For a third-order accurate three-step scheme, we have indices 
i = 1,2,3. In order to advance from oj h at time t n to of +l 



In this section, the schemes under consideration are vali- 
dated on a square domain by solving the Taylor-Green decaying 
vortex problem. The Taylor-Green vortex problem is the two- 
dimensional, unsteady flow of a decaying vortex, which is an 
exact closed form solution to the incompressible Navier-Stokes 
equations in Cartesian coordinates. The exact solution of this 
vortex flow in a [0, 2n] x [0, 2n] domain with periodic boundary 
conditions is given by 



oj e (x,y,t) = 2KCOs(Kx)cos(Ky)e 2l<2t/Re 



(42) 



where k is an integer which represents the number of vortices 
in each direction. In order to quantify the effective order of 
accuracy for each scheme, using the difference between exact 
and computed solutions, we compute the discrete L 2 norm as 



\M\l 2 



1 



Af, Ny 



(43) 



In order to test the spatial convergence rate, we first perform 
numerical simulations for Re = 1 and k = 4 with all the spatial 
schemes introduced earlier using the third-order TVD Runge- 
Kutta scheme with At = 10" 4 . We choose this small time step 
(wherein the maximum CFL number is around 0.002 for the 
128 2 resolution case) to eliminate possible temporal discretiza- 
tion errors due to the Runge-Kutta time integration schemes. A 
study with At = 2 x 10" 4 , 10" 4 , and 5 x 10" 5 shows that the pre- 
dicted flow field is independent of At. To further verify whether 
the simulations are independent of the errors associated with 
temporal discretization, we also run the same computations us- 
ing the fourth-order Runge Kutta scheme, which yields exactly 
the same results. 
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Table 2: Discrete L2 error norms and corresponding convergence rates showing the effective order of the spatial difference schemes for the Taylor-Green vortex 
problem at t = 0.1 for Re = 1 and k = 4, obtained using the time step At = 10 -4 . The final column shows the CPU times in seconds for the 128 2 resolution 
computations. 



Scheme 


IMI L2 


n 32 2 


H< 


32 2 

n 


n< 2 


n 641 
n m 2 


IMlif 


CPU 


pseudospectral method (PS) 


0.76E-9 




6.63E-9 




1.61E-8 




2.36E-8 


210 


6th-order explicit difference (ED6) 


1.04E-2 


5.67 


2.05E-4 


5.92 


3.39E-6 


6.00 


5.28E-8 


46 


4th-order explicit difference (ED4) 


3.28E-2 


3.92 


2.16E-3 


3.97 


1.38E-4 


4.00 


8.65E-6 


44 


2nd-order explicit difference (ED2) 


1.44E-1 


2.29 


2.94E-2 


2.09 


6.91E-3 


2.03 


1.70E-3 


42 


6th-order compact difference (CD6) 


2.90E-3 


6.19 


3.98E-5 


6.06 


5.96E-7 


6.14 


8.47E-9 


79 


4th-order compact difference (CD4) 


1.58E-2 


4.17 


8.74E-4 


4.05 


5.28E-5 


4.02 


3.26E-6 


76 


4th-order Arakawa scheme (A4) 


3.28E-2 


3.92 


2.16E-3 


3.97 


1.38E-4 


4.00 


8.65E-6 


52 


2nd-order Arakawa scheme (A2) 


1.44E-1 


2.29 


2.94E-2 


2.09 


6.91E-3 


2.03 


1.70E-3 


45 


4th-order DRP scheme (DRP4) 


3.28E-2 


3.92 


2.16E-3 


3.97 


1.38E-4 


4.00 


8.65E-6 


45 
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Figure 2: Convergence of the discrete L2 error norms at time t = 0.1 for the 
Taylor-Green vortex decaying problem at Re = 1 . 



In evaluating the performance of different finite difference 
formulations, the computed discrete L2 norms at time t = 0.1 
(when a 95% decay of the initial field has happened) are tab- 
ulated in Tableland are also plotted in Fig. [2] Results show 
that the theoretical order of accuracies for all the schemes are 
obtained in practice. A linear reduction rate is obtained with 
the slope of all curves being their theoretical values, verifying 
the correct spatial convergence rate for each scheme. The most 
accurate difference scheme is found to be the sixth-order com- 
pact difference scheme. We also demonstrate that the errors 
associated with all the high-order formulations are several or- 
der magnitudes lower than the second-order formulations. It 
can be seen from Fig. [2] that the predicted result of the 6th- 
order compact scheme becomes more accurate than the pseu- 
dospectral solution for the case with 128 2 resolution. This is 
mainly due to the accumulation of round off errors in the FFTs 



of the pseudospectral method. It can also be seen that the as- 
sociated error increases slightly with increasing resolution for 
pseudospectral computations, again, due to the round-off errors 
associated with the FFTs. That is actually why Lele [ 35 ] called 
these compact schemes "spectral-like schemes." It can also be 
seen that the compact schemes give more accurate results than 
the explicit finite difference schemes due to their smaller trunca- 
tion error. Since the viscous difference formulations are equiv- 
alent, it should be noted here that the fourth-order Arakawa, 
dispersion-relation-preserving, and explicit difference schemes 
predict almost the same answer because of the inherent dynam- 
ics of the Taylor-Green vortex validation problem. This prob- 
lem does not have a convective flow field, hence, the center of 
the each vortex stays stationary. 

Next, we perform a convergence rate analysis for the tempo- 
ral schemes. As we can see from the previous analysis, the 
pseudospectral scheme provides much more accurate results 
than the finite difference schemes, especially for lower reso- 
lutions. In order to be able to test the convergence rate of the 
temporal schemes properly, we choose our spatial discretization 
scheme to be the pseudospectral scheme with a resolution of 
16 2 , for which the temporal discretization error becomes greater 
than the spatial discretization error. The computed error norms 
are shown in Table |3] at time t = 20 for Re = 1000. It can be 
seen from Fig.[3]that a linear reduction rate is obtained with the 
slope of all curves being 3 for all the third-order Runge-Kutta 
schemes, and 4 for the fourth-order Runge-Kutta scheme. We 
report that all the third-order Runge-Kutta schemes utilized in 
this study produce the same result within negligible differences 
due to round-off error. 

The Taylor-Green vortex problem is one of the simplest in- 
compressible flow test problems, having a smoothly decaying 
field. We solve it to be able to perform an exact error conver- 
gence rate analysis by comparing our results to the analytical 
solution of the problem. In the following sections, the effects 
of the order of the schemes and their global conservation prop- 
erties will be analyzed for more challenging convective dom- 
inated flow problems. A truly fair comparison might require 
implementation of the same discretization scheme for the vis- 
cous terms. Therefore, in order to do a fair comparison among 
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Table 3: Discrete L2 error norms and corresponding convergence rates showing the effective order of the temporal integration schemes for the Taylor-Green vortex 
problem at t = 20 for Re = 1000 and k = 4, obtained using the pseudospectral scheme with a spatial resolution of 16 2 . 



Scheme IMIf =ai rate IMIf =a2 rate |Mlf =a4 rate |M 



4th-order Runge-Kutta (RK4) 


1.23E-12 


4.03 


2.01E-11 


4.01 


3.23E-10 


4.02 


5.23E-09 


TVD Runge-Kutta (TVDRK3) 


1.96E-09 


3.00 


1.57E-08 


3.01 


1.26E-07 


3.01 


1.02E-06 


Symmetric (SYMRK3) 


1.96E-09 


3.00 


1.57E-08 


3.01 


1.26E-07 


3.01 


1.02E-06 


Predictor/Corrector (P/CRK3) 


1.96E-09 


3.00 


1.57E-08 


3.01 


1.26E-07 


3.01 


1.02E-06 


Inhomogeneous (INHRK3) 


1.96E-09 


3.00 


1.57E-08 


3.01 


1.26E-07 


3.01 


1.02E-06 


Williamson (WILRK3) 


1.96E-09 


3.00 


1.57E-08 


3.01 


1.26E-07 


3.01 


1.02E-06 



10" 5 r ■ RK4 
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where the constants cr and 6 determine the thickness of the shear 
layer and the amplitude of the initial perturbation, respectively. 
In our computations, the thickness parameter is cr = 15 /n, the 
perturbation amplitude is S = 0.05, the Reynolds number is 
Re = 10 4 , and the time step is At = 10" 3 . We choose this 
time step (wherein the corresponding maximum CFL number 
is 0.25 for the case with a resolution of 1024 2 ) to concentrate 
on the differences among the spatial discretization schemes by 
eliminating the temporal integration errors in the simulations. 
Since the dominant error is due to the spatial difference opera- 
tor rather than the temporal integration method, we employ the 
TVD Runge-Kutta scheme for the time advancement. We also 
perform, but do not report in detail here, a time step refinement 
study revealing that the temporal discretization error is indeed 
negligible. Moreover, we use the same discrete approximation, 
the 6th-order compact difference scheme, for the viscous terms 
in all the cases. In the following, we investigate the behavior of 
the spatial discretization schemes for convective terms by solv- 
ing the double shear layer problem for various resolutions. 



Figure 3: Convergence of the discrete L2 error norms for the temporal schemes 
at time t = 20 for the Taylor-Green vortex decaying vortex problem at Re = 
1000. 



spatial difference schemes, the viscous diffusion terms are ap- 
proximated using the six-order accurate compact scheme for all 
the finite difference schemes being considered for the rest of the 
paper. 

5. Double shear layer problem 

The double shear layer problem was introduced by Bell, 
Colella and Glaz l52l to test the performance of a projection 
method, and was studied later by Minion and Brown l53l us- 
ing various schemes to determine the effects of the grid resolu- 
tion on solutions of the unsteady, incompressible Navier-Stokes 
equations. It is a benchmark problem for testing the accuracy 
and resolution of a time dependent numerical method. The dou- 
ble shear layer problem with periodic boundary conditions in a 
square domain [0, 2n] x [0, 2n] is subjected to the following ini- 
tial conditions 




oj(x, y, 0) 



6 cos(x) - cr cosh 2 [cr(y - f )] if y < 71 
6 cosO) + cr cosh" 2 [<T(Tf - y)] if y > n 



Figure 5: The three-dimensional view of the shear layers at time t - 10 by 
plotting the magnitude of the vorticity as the z-axis. 








Figure 4: Evolution of the double shear layer problem for Re - 
domain. 



10 with the perturbation and thickness parameters of 5 = 0.05 and cr = 15 /n in a doubly periodic 



Table 4: Discrete L2 error norms and computational efficiencies of the spatial schemes for the double shear layer problem at t ■ 
solution for computing the L2 norm is obtained using the pseudospectral method with a resolution of 1024 2 . 



10 for Re = 10 . The reference 



Scheme 


IMgf 


CPU 2562 (hrs) 


iMiif 


CPU 5122 (hrs) 


IMC 4 ' 


CPlW(hrs) 


PS 


5.91E-3 


2.94 


3.34E-6 


17.65 




97.03 


ED6 


8.10E-2 


0.67 


6.86E-3 


3.08 


1.87E-4 


19.29 


ED4 


1.48E-1 


0.66 


2.53E-2 


2.99 


2.04E-3 


17.84 


ED2 


3.97E-1 


0.64 


1.81E-1 


2.95 


5.86E-2 


16.93 


CD6 


3.36E-2 


0.92 


1.03E-3 


4.10 


8.10E-5 


23.60 


CD4 


6.93E-2 


0.90 


6.47E-3 


4.05 


3.71E-4 


21.73 


A4 


1.45E-1 


0.76 


2.52E-2 


3.37 


2.04E-3 


19.32 


A2 


4.45E-1 


0.69 


2.01E-1 


3.15 


6.46E-2 


17.92 


DRP4 


5.69E-2 


0.67 


8.23E-3 


3.05 


1.02E-3 


18.07 



The perturbed shear layers roll up into a single vortex and 
the shear layers become thinner and thinner as the flow evolves 
over time, as shown in Fig. [4] The three-dimensional view of 
the shear layers developed at time t = 10 is also illustrated in 
Fig. |5j which shows the magnitude of the vorticity as the z- 
axis. This shear layer problem is particularly important in that 



the presence of thinner and thinner shear layers as the flow field 
evolves in time is not captured by low grid resolution repre- 
sentations. The Gibbs phenomenon, numerical oscillations oc- 
curring near sharp vorticity gradients, occurs for underresolved 
simulations in which the grid size is larger than the shear layer 
thickness. Since sharp discontinuities are observed in the vor- 
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ticity field, this benchmark flow problem is very appropriate to 
test the overall accuracy of the methods presented in this study. 

Vorticity contours at t = 10 obtained using all the spatial dis- 
cretization schemes introduced earlier are shown in Figs. [6|8] 
for different spatial resolutions. To make the distinctions among 
the schemes more clear, the centerline vorticity distributions 
along the v-axis at t = 10 are also shown in Fig. 9|[TT for differ- 
ent resolutions. Our comparison in these figures includes four 
families of fourth-order difference schemes, two second-order 
schemes, and two sixth-order difference schemes, as well as the 
pseudospectral scheme. Fig. [6] illustrates a comparison of vor- 
ticity fields obtained by using 25 6 2 resolution for these nine dif- 
ferent spatial discretization schemes. The corresponding com- 
parison of centerline vorticity distributions is also provided in 
Fig. [9] As we can see from these figures, both the sixth-order 
compact scheme and the pseudospectral scheme yield almost 
the same accuracy. It can also be seen that all the high-order for- 
mulations provide much better accuracy than the second-order 
schemes. A truly fair comparison among all the finite difference 
formulations presented in this study is to look at the fourth- 
order discretization schemes. Results show that there is no 
significant difference between explicit difference and Arakawa 
schemes. The discrete conservation properties of the Arakawa 
scheme do not improve the overall accuracy of the whole so- 
lution procedure in this case. On the other hand, the compact 
difference and dispersion-relation-preserving schemes result in 
better accuracy than the Arakawa or explicit difference schemes 
by eliminating numerical oscillations. Increasing the resolu- 
tion to 51 2 2 , as illustrated in Fig. [7]and Fig. 10 we show that 
it is possible to obtain numerical oscillation-free results only 
by using high-order schemes. Second-order schemes, both the 
explicit difference and the Arakawa scheme, still produce nu- 
merical oscillations. Finally, we perform computations by us- 
ing 1024 2 resolution and results are demonstrated in Fig. [8] and 
Fig. 1 1 for vorticity contour fields and centerline vorticity dis- 
tributions, respectively. We can see that all numerical schemes 
provide an accurate oscillation-free solution to the double shear 
layer problem with this resolution. 

As a result of this study, a rule of thumb is that second-order 
central difference schemes, whether standard explicit difference 
or conservative Arakawa scheme, require about twice resolu- 
tion (in each direction) to achieve the same results obtained by 
a high-order finite difference or pseudospectral scheme. There- 
fore, this study demonstrates that the numerical oscillations 
that occur in underresolved simulations with strong shear lay- 
ers can be eliminated efficiently by using the high-order formu- 
lations with the same resolution. For example, at a relatively 
high Reynolds number Re = 10 4 , it can clearly be seen that 
a resolution of 512 2 is not high enough to capture the correct 
physics without numerical oscillations using a low-order finite 
difference formulation. To be able to obtain all the small-scale 
physical layers, at least a 1024 2 resolution needs to be used 
with a spatially second-order accurate scheme for solving un- 
derlying equations. On the other hand, we demonstrated that 
the properly-resolved physical behavior without a Gibbs phe- 
nomenon is obtained using a resolution of 512 2 with only a 
slight increase in the computational price. 



In evaluating the efficiency of the different formulations, the 
CPU times are tabulated in Table [4] In this table, we also 
present the corresponding discrete L2 norms which show the 
average root-mean-square deviations of the vorticity field from 
the data obtained by well-resolved pseudospectral calculation 
using 1024 2 resolution. These results demonstrate that high fi- 
delity numerical simulations can be obtained using the high- 
order finite difference methods with a speed-up factor of 5 if 
we compare them with the pseudospectral computation. More 
importantly, if we compare them with their second-order coun- 
terparts, the additional cost due to high-order formulations be- 
comes quite small. The data presented in Table [4] also demon- 
strates that the computational speed-up increases with increas- 
ing spatial resolution as expected. With this in mind, one of 
the interesting results is that the errors associated with both 
the second-order explicit difference (ED2) and the second-order 
Arakawa (A2) schemes with 1024 2 resolution are greater than 
the error associated with the sixth-order compact scheme (CD6) 
with a much smaller 25 6 2 resolution. This clearly demonstrates 
that the high-order difference formulations can provide an ef- 
ficient way of solving flow problems even with strong shear 
components. 
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Figure 9: The centerline vorticity distributions at x = n for the double shear 
layer problem, plotted for a resolution of 256 2 . Comparison of the different 
numerical schemes at time t = 10. Inset: close-up of boxed area. 



6. Two-dimensional decaying turbulence 

Two-dimensional homogenous decaying turbulence is an in- 
compressible flow problem in which the kinetic energy decays 
l54l . The problem of the free decay of two-dimensional tur- 
bulence is to determine how abundant populations of vortices 
freely evolve with time JT4l . In this study, we solve this prob- 
lem to examine the accuracy and efficiency of the presented fi- 
nite difference approximations. The results are compared with 
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Figure 6: Comparison of the numerical schemes for the double shear layer problem at time t = 10 with a resolution of 256 2 . (a) pseudospectral (PS) method, (b) 
sixth-order explicit difference (ED6) method, (c) fourth-order explicit difference (ED4) method, (d) second-order explicit difference (ED2) method, (e) sixth-order 
compact difference (CD6) method, (f) fourth-order compact difference (CD4) method, (g) fourth-order Arakawa (A4) method, (h) second-order Arakawa (A2) 
method, and (i) fourth-order dispersion-relation-preserving (DRP4) method. The vorticity contour layouts are identical in all nine cases illustrating 27 equidistant 
levels in the interval [-4.5, 4.5]. 
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Figure 7: Comparison of the numerical schemes for the double shear layer problem at time t = 10 with a resolution of 512 2 . (a) pseudospectral (PS) method, (b) 
sixth-order explicit difference (ED6) method, (c) fourth-order explicit difference (ED4) method, (d) second-order explicit difference (ED2) method, (e) sixth-order 
compact difference (CD6) method, (f) fourth-order compact difference (CD4) method, (g) fourth-order Arakawa (A4) method, (h) second-order Arakawa (A2) 
method, and (i) fourth-order dispersion-relation-preserving (DRP4) method. The vorticity contour layouts are identical in all nine cases illustrating 27 equidistant 
levels in the interval [-4.5, 4.5]. 
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Figure 8: Comparison of the numerical schemes for the double shear layer problem at time t = 10 with a resolution of 1024 2 . (a) pseudospectral (PS) method, (b) 
sixth-order explicit difference (ED6) method, (c) fourth-order explicit difference (ED4) method, (d) second-order explicit difference (ED2) method, (e) sixth-order 
compact difference (CD6) method, (f) fourth-order compact difference (CD4) method, (g) fourth-order Arakawa (A4) method, (h) second-order Arakawa (A2) 
method, and (i) fourth-order dispersion-relation-preserving (DRP4) method. The vorticity contour layouts are identical in all nine cases illustrating 27 equidistant 
levels in the interval [-4.5, 4.5]. 
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Figure 10: The centerline vorticity distributions at x = n for the double shear 
layer problem, plotted for a resolution of 512 2 . Comparison of the different 
numerical schemes at time t = 10. Inset: close-up of boxed area. 




Figure 1 1 : The centerline vorticity distributions at x = n for the double shear 
layer problem, plotted for a resolution of 1024 2 . Comparison of the different 
numerical schemes at time t = 10. Inset: close-up of boxed area. 



those of the pseudospectral method. According to the KBL the- 
ory for two-dimensional turbulence this system has an inertial 
range in the energy spectrum that is proportional to k~ 3 in the 
inviscid limit. 

The computational domain is a square with sides of length 
2n. Periodic boundary conditions are applied. The initial en- 




Log (k) 

Figure 12: Initial energy spectrum for the two-dimensional decaying turbulence 
problem. 



ergy spectrum in Fourier space is given by 

^ 2s+l 

exp 



2 \k, 



(45) 



where k = |k| = + ky. The maximum value of initial en- 
ergy spectrum occurs at wavenumber k p which is assumed to be 
k p = 12 in this study. The coefficient a s normalizes the initial 
kinetic energy and is given by 



a s = 



(2s + 1) 

2 s si 



s+l 



(46) 



where s is a shape parameter. In this study, we take s = 3. The 



corresponding energy spectrum given in Eq. ( 45 ) is illustrated 
in Fig. [12] The magnitude of vorticity Fourier coefficients re- 
lated to the assumed initial energy spectrum becomes 



-E(k) 



(47) 



The initial vorticity distribution in Fourier space is then ob- 
tained by introducing a random phase 



-E(k) ^ (k) 

71 



(48) 



where the phase function is given by £(k) = £(k) + rj(k), where 
£(k) and r](k) are independent random values chosen in [0, 2tt] 
at each coordinate point in the first quadrant of the k x -k y plane. 
The conjugate relations for other quadrants are 

^(~k x > ky) = —^(k x , ky) 

K-k X9 -ky) = -€(k X9 ky) 
{(k X ,-ky)={(k X ,ky) 

rj(-k x , k y ) = T](k x , k y ) 
Tj(-k x ,-ky) = -rj(k x ,k y ) 
rj(k x ,-ky) = -r 1 (k x ,ky). (49) 

After the randomization process described above for the 
phases, the initial vorticity distribution in the physical space is 
computed by taking the inverse FFT. However, we use a random 
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phase generator to obtain initial flow field, it should be noted 
that we will use exactly same initial vorticity field when we 
compare the numerical schemes for the same resolution. The 
corresponding initial vorticity distribution is a vortex popula- 
tion which satisfies the divergence free condition and provides 



the energy spectrum given in Eq. (45 ). In the following text, we 
first investigate how homogeneous two-dimensional decaying 
turbulence varies with Re when using the sixth-order compact 
difference scheme. Then we will focus our attention on the be- 
havior of a selection of numerical methods for computations of 
two-dimensional turbulence. It should be noted here that the 
diffusion terms are approximated using the sixth-order compact 
scheme for all difference formulations in order to provide a fair 
comparison between them. To make the simulation independent 
of the errors associated with temporal discretization, a time step 
of At = 2 x 10" 4 is used in all simulations except when the tem- 
poral discretization schemes are evaluated by a time refinement 
analysis. A study with At = 1 x 10" 4 , 2 x 10" 4 , and 4 x 10" 4 
also show that the predicted flow field is independent of At. 

6.1. Reynolds number dependence 

Starting from a divergence free initial vorticity field having 
the specified energy spectrum, we solve Eq. ([l} for various 
Reynolds numbers. Fig. 13 shows the evolution of the vortic- 
ity field for Re = 1000. As shown in Fig. [T3] initial flow field 
at time t = consists of many vortices according to the spec- 



trum given by Eq. (45 ). As we can see, a filamentation process 
occurs, initially randomly distributed vortices start to interact 
with each other and become bigger vortices with time by a vor- 
tex merging mechanism. As also illustrated in Fig. [13] as time 



evolves, the exterior strain field from each vortex deforms the 
vorticity field of the other one so that the vorticity fields wrap 
around each other. 

In order to examine the characteristics of two-dimensional 
decaying turbulence simulations and their dependencies on Re, 
we first define two statistical measures, one is the energy spec- 
trum in wave space, and the other is the structure function in 
physical space. The energy spectrum is defined as 



E(Kt)=^k 2 \HKt)\ 2 

and the angle averaged energy spectrum is 

E(k, t)= Yj 

k<\k\<k+i 



(50) 



(51) 



It is known from the KBL theory that the energy spectrum in the 
inertial range approaches the classical k~ 3 scaling in the limit 
of infinite Reynolds number. On the other hand, the statistics of 
two-dimensional turbulent flow can be further investigated con- 
sidering powers of vorticity differences in the physical space. 
A commonly used statistical quantity in two-dimensional tur- 
bulence is the second-order vorticity structure function which 
is defined as 



with r = |r| being the spatial separation. Statistical properties 
of decaying two-dimensional turbulence are investigated here 
by numerical simulations for different Reynolds numbers. 

In this part of our analysis, in order to elucidate the flow field 
evolution for various Re and how the energy spectra and struc- 
ture functions depend on Re, we perform a set of numerical 
experiments for Re = 100, 300, 600, 1000, 3000, and 6000. 
All numerical simulations are performed using a resolution of 
1024 2 and start from the same initial condition. Fig. 
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demon- 
strates the evolution of the angle averaged energy spectrum for 
different Reynolds numbers. As we can see, the energy spec- 
trum becomes steeper than k~ 3 for Re = 100 due to the smaller 
filamentation and less interaction among vortices. However, 
increasing Re, the energy spectrum in the inertial range ap- 
proaches towards the classical k~ 3 limit, in agreement with KBL 
theory of two-dimensional turbulence. Similarly, the evolution 
of the second-order structure functions for various Reynolds 

showing r 2 scaling for small 



15 



numbers is exhibited in Fig. 
separations in all the cases. They flatten for large separations 
as predicted by the KBL theory of two-dimensional turbulence 
in the inviscid limit. Comparisons of the angle averaged en- 
ergy spectra and the second-order vorticity structure functions 



at time t = 2 are also shown in Fig. [FT] and Fig. [18] respec- 
tively. As demonstrated by Fig. [17] the angle averaged energy 
spectrum defined by Eq. (51) asymptotically reaches the k~ 3 
spectrum in the inertial range as Re increases. We find that the 
Reynolds number dependency is more stringent if we look at 
the turbulence statistics in wave space such the angle averaged 
energy spectrum. The corresponding instantaneous vorticity 



fields at time t = 2 are compared in Fig. 16 for the same set 



of Reynolds numbers. As we can see from Fig. 16 the amount 
of filamentation increases for higher Reynolds numbers. Due to 
the smaller convection in lower Reynolds numbers, the interac- 
tion between two vortices is not as strong as that of the higher 
Reynolds number cases, such that less amount of filamentation 
occurs and results in bigger vortices with less magnitude of vor- 
ticity during the evolution because of faster decay rate. 

Noticeably, our results provide numerical proof of the 
KBL theory for decaying two-dimensional turbulence using 
the sixth-order compact finite difference scheme for spatial 
discretization along with the third-order TVD Runge-Kutta 
scheme for the time advancement. However, we emphasize that 
our primary goal in this study is to investigate the behavior of 
different numerical schemes for long time integration of turbu- 
lence simulations. In the following, we first test the temporal 
schemes and their numerical stabilities, and then focus on the 
evaluation of the behavior of spatial discretization schemes for 
computations of decaying turbulence. 

6.2. Comparison of temporal schemes 

We devote this section to investigating the characteristics of 



the temporal integration schemes described in Section 3.7 In 



«r)V(Wx + r)-^(x)| 2 ) 



(52) 



general, the errors associated with temporal discretization are 
smaller than the errors associated with the spatial discretiza- 
tion, especially if we use high-order methods. Although time 
step limitations in an explicit time integration method can be an 
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Figure 13: Evolution of the vorticity field in decaying turbulence for Re = 1000. Initially randomly distributed vortices start to interact with each other and merge 
to form larger vortices with time. 
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Figure 14: Evolution of energy spectra in decaying turbulence for (a) Re=100, (b) Re=300, (c) Re=600, (d) /te=1000, (e) tfe=3000, and (f) /te=6000. The energy 
spectrum in the inertial range flattens towards the classical k~ 3 scaling limit as Re increases, in agreement with the KBL theory of two-dimensional turbulence. 
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Figure 15: Evolution of the second-order vorticity structure functions in decaying turbulence for (a) 7te=100, (b) 7te=300, (c) /te=600, (d) 7te=1000, (e) 7te=3000, 
and (f) /te=6000. The line r 2 is included in each subfigure. 




Figure 17: Reynolds number dependence of the energy spectrum in wave space. 
In the inertial range, the angle averaged energy spectrum, defined by Eq. {5l) , 
asymptotically reaches k~ 3 scaling with increasing Reynolds number. 



Figure 18: Reynolds number dependence of the second-order vorticity struc- 
ture functions in physical space showing the second-order vorticity structure 
function, as defined by Eq. {52) . 



impediment to efficient solutions, for turbulent flow computa- 
tions small time steps are indeed desirable to resolve the small 
turbulent time scales. In fact, a balanced spatio-temporal dis- 
cretization dictates a CFL number of 0(1) and large time steps 
can have the same undesirable effect on the solution as a coarse 
spatial resolution or a low-order dissipative spatial discretiza- 
tion. 

In this section, we evaluate the behavior of different explicit 
Runge-Kutta schemes introduced earlier for time integration 
procedure by performing a time step refinement study. The do- 
main of absolute stability for the explicit Runge-Kutta schemes 
for scalar equations can be found in many references including 
Canuto et al. l55l and Buthcher [56]. What is true for scalar 
equations is not necessarily true for systems as advocated by 



Pruett et al. |47l . Our method of estimating stability limits is 
somewhat heuristic here, but it becomes a useful guideline in 
practical applications for a wide variety of numerical methods. 

For this purpose, we first compare the accuracy and stability 
of the temporal schemes by performing a time refinement study 
using the the sixth-order compact difference scheme as a spa- 
tial discretization method with 512 2 resolution. Table [U shows 
the discrete L2 error norms at t = 10 for Re = 10 3 obtained 
by the explicit Runge-Kutta schemes introduced earlier varying 
the time steps. The corresponding maximum CFL numbers for 
each computation are also documented in the table. Deviations 
of oj are computed from the base solution which is obtained by 
the forth-order Runge-Kutta method with a small time step of 
At = 2 x 10" 4 . In terms of numerical stability, it is demon- 
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Figure 16: Vorticity contour plots at time t = 2 for various Reynolds numbers. Stronger filamentation processes occur with increasing Reynolds number. 



strated in the table that the fourth-order Runge-Kutta scheme 
has a greater stability region than the third-order ones as we 
expect. Among the third-order Runge-Kutta schemes for the 
explicit time integration, the optimal TVD Runge-Kutta scheme 
becomes more accurate than all the low-storage schemes. How- 
ever, the difference between the Williamson scheme and the 
TVD Runge-Kutta scheme is negligible. On the other hand, we 
find that the predictor/corrector type Runge-Kutta scheme ex- 
hibit substantially lower accuracy than other third-order meth- 
ods. 

Next, we perform a similar analysis by using the second- 
order conservative Arakawa schemes for the same parameters. 
Results are demonstrated in Table [6] showing the similar trends 
as observed in the Table [5] for the compact scheme. The in- 
teresting outcome is that the second-order scheme can run us- 
ing almost triple the effective time step if we compare with 
the sixth-order scheme due to a wider stability region. Al- 
though the results are not shown here, we also performed sim- 
ilar studies for the fourth-order Arakawa and the fourth-order 
compact schemes. Between the second-order scheme and the 



fourth-order ones, the second-order spatial operator runs at 
nearly double the effective time step. These results suggest 
that the second-order scheme becomes more competitive with 
other high-order accurate spatial schemes by allowing larger 
time steps. On the other hand, as we will show in the follow- 
ing section, the second-order scheme usually requires double 
the resolution in each direction to be able to obtain a similar 
accuracy. Even if using a low order scheme allows us to use 
larger time steps, still, the high-order schemes become more 
effective in terms of the tradeoff between the accuracy and ef- 
ficiency. We also demonstrate that although we use an optimal 
fast Poisson solver, the CPU time increases approximately 8 
times by doubling the resolution in each direction. In fact, it 
would be larger if we use other suboptimal Poisson solver. In 
other worlds, that fraction of CPU time is much bigger in other 
geometries in which the solution of the elliptic equation domi- 
nates the computational time. It should also be noted that the ef- 
ficiency of the high-order schemes would be more pronounced 
for three dimensional computations. 
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Table 5: Discrete L2 error norms for the decaying turbulence problem at t = 10 for Re = 10 3 using the sixth-order compact difference scheme as a spatial 
discretization method. The reference solution for computing the L2 norm is obtained using the RK4 method with a time step of At = 2 x 10 -4 . A dash indicates the 
numerical instability. 

Scheme IMI^ 1 IMI™' 4 IM|™« IM|™-° IMI^™ 5 

(A? = 4 x 10" 4 ) (Ar = 8 x 10" 4 ) (At = 16 x 10" 4 ) (At = 32 x 10" 4 ) (At = 40 x 10" 4 ) (Af = 50 x 10" 4 ) 



RK4 


7.50E-7 


4.01E-6 


5.31E-5 


7.87E-4 


1.84E-3 




TVDRK3 


3.92E-5 


3.11E-4 


2.46E-3 


1.94E-2 






SYMRK3 


4.65E-4 


1.94E-3 


8.46E-3 


3.98E-2 






P/CRK3 


1.32E-3 


5.23E-3 


2.06E-2 


8.16E-2 






INHRK3 


4.02E-5 


3.18E-4 


2.51E-3 


1.96E-2 






WILRK3 


3.99E-5 


3.14E-4 


2.49E-3 


1.95E-2 







Table 6: Discrete L2 error norms for the decaying turbulence problem at t = 10 for Re = 10 3 using the second-order conservative Arakawa scheme as a spatial 
discretization method. The reference solution for computing the L2 norm is obtained using the RK4 method with a time step of At = 2 x 10 -4 . A dash indicates the 
numerical instability. 



Scheme IMI™' 2 INIg™ N^f™ IMI^" IM^™ IMIg™ 





(At = 8 x 1(T 4 ) 


(At = 16 x 10- 4 ) 


(At = 50 x 1(T 4 ) 


(At = 100 x 10" 4 ) 


(At = 125 x 10~ 4 ) (At = 140 x 10~ 4 ) 


RK4 


2.05E-6 


1.94E-5 


1.76E-3 


3.31E-2 


9.89E-2 


TVDRK3 


1.65E-4 


1.33E-3 


4.34E-2 


4.16E-1 




SYMRK3 


2.05E-3 


8.61E-3 


1.04E-1 


5.70E-1 




P/CRK3 


5.78E-3 


2.28E-2 


2.09E-1 


7.18E-1 




INHRK3 


1.77E-4 


1.41E-3 


4.52E-2 


4.20E-1 




WILRK3 


1.75E-4 


1.40E-3 


4.47E-2 


4.16E-1 





6.3. Comparison of spatial schemes 

Our main objective is to investigate and analyze the behaviors 
of the spatial numerical schemes presented in Section|3]for two- 
dimensional decaying turbulence simulations. In this section, 
four different families of high-order accurate finite difference 
representations are compared with pseudo spectral simulations 
for various Reynolds numbers. We also compare high-order ac- 
curate schemes with the second-order accurate schemes. From 
our previous time refinement study, we conclude that the trun- 
cation error produced by Runge-Kutta time integration schemes 
is much smaller than the discretization schemes when we use a 
time step within the stability region. The TVD-RK scheme with 
At = 2 x 10" 4 is used as a time integration algorithm for the 
following analysis in which the CFL number becomes O(0. 1). 
Therefore, we can compare different spatial schemes without 
contaminating the solution due to time stepping algorithm. We 
also highlight that we use the sixth-order compact difference 
formulation for discretization of viscous terms in all schemes 
tested here. 

First, we present a set of numerical experiments using a reso- 
lution of 1024 2 for Re = 1000. The vorticity fields obtained by 
nine different spatial schemes at time t = 5 are presented and 



5r 



compared in Fig. 19) A comparison of the vorticity distributions 
on the vertical centerline is also illustrated in Fig. [23] by show- 
ing also a close-up figure to make the comparison more clear. 
Furthermore, a comparison of statistical behavior in terms of 




Figure 23: Comparison of the numerical schemes for Re = 1000 with a resolu- 
tion of 1024 2 (Re c = 6.13) at time t = 5. The centerline vorticity distributions 
at x = n are plotted. Inset: close-up of boxed area. 



the second-order vorticity structure functions is also demon- 
As we can see, the high-order schemes 



strated in Fig. 24 
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Figure 19: Comparison of the numerical schemes at time t = 5 for Re = 1000 with a resolution of 1024 2 (Re c = 6. 13). (a) pseudospectral(PS) method, (b) sixth-order 
explicit difference (ED6) method, (c) fourth-order explicit difference (ED4) method, (d) second-order explicit difference (ED2) method, (e) sixth-order compact 
difference (CD6) method, (f) fourth-order compact difference (CD4) method, (g) fourth-order Arakawa (A4) method, (h) second-order Arakawa (A2) method, and 
(i) fourth-order dispersion-relation-preserving (DRP4) method. The vorticity contour layouts are identical in all nine cases illustrating 27 equidistant levels in the 
interval [-13, 13]. 
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Figure 20: Comparison of the numerical schemes at time t = 10 for Re = 1000 with a resolution of 512 2 (Re c = 12.27). (a) pseudospectral(PS) method, (b) 
sixth-order explicit difference (ED6) method, (c) fourth-order explicit difference (ED4) method, (d) second-order explicit difference (ED2) method, (e) sixth-order 
compact difference (CD6) method, (f) fourth-order compact difference (CD4) method, (g) fourth-order Arakawa (A4) method, (h) second-order Arakawa (A2) 
method, and (i) fourth-order dispersion-relation-preserving (DRP4) method. The vorticity contour layouts are identical in all nine cases illustrating 41 equidistant 
levels in the interval [-8, 8]. 
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Figure 21 : Comparison of the numerical schemes at time t = 6 for Re = 1000 with a resolution of 256 2 (Re c = 24.54). (a) pseudospectral(PS) method, (b) sixth-order 
explicit difference (ED6) method, (c) fourth-order explicit difference (ED4) method, (d) second-order explicit difference (ED2) method, (e) sixth-order compact 
difference (CD6) method, (f) fourth-order compact difference (CD4) method, (g) fourth-order Arakawa (A4) method, (h) second-order Arakawa (A2) method, and 
(i) fourth-order dispersion-relation-preserving (DRP4) method. The vorticity contour layouts are identical in all nine cases illustrating 23 equidistant levels in the 
interval [-11, 11]. 
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Figure 22: Comparison of the numerical schemes at time t = 10 for Re = 100 with a resolution of 512 2 (Re c = 1.23). (a) pseudospectral(PS) method, (b) sixth-order 
explicit difference (ED6) method, (c) fourth-order explicit difference (ED4) method, (d) second-order explicit difference (ED2) method, (e) sixth-order compact 
difference (CD6) method, (f) fourth-order compact difference (CD4) method, (g) fourth-order Arakawa (A4) method, (h) second-order Arakawa (A2) method, and 
(i) fourth-order dispersion-relation-preserving (DRP4) method. The vorticity contour layouts are identical in all nine cases illustrating 19 equidistant levels in the 
interval [-0.45,0.45]. 
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and the pseudo spectral schemes yield similar flow fields. The 
second-order accurate schemes, both the explicit difference and 
the Arakawa schemes, predict slightly less accurate results. The 
sixth-order compact scheme exhibits the best prediction capa- 
bility among the different formulations tested. In evaluating the 
performance of the different finite difference approximations, 
the deviation of vorticity field from the spectral solution is cal- 
culated by using the definition of discrete L2 norm given by 
Eq. ( |49| ). Table [7] shows the L 2 norms, CPU times, and corre- 
sponding speed-up ratios for the finite difference formulations. 
The speed-up ratio is defined here as the ratio of CPU times for 
the computations performed using the pseudospectral method 
and those performed using finite difference schemes. First, all 
the finite difference methods are considerably more efficient 
than the pseudospectral method. The data in the table exhibits 
that the high-order schemes provide the same accuracy level 
as the Fourier-Galerkin pseudospectral method, while signifi- 
cantly decreasing its computational cost with a speed-up factor 
of 5 over the pseudospectral method. Among the fourth-order 
schemes, both the compact and dispersion-relation-preserving 
schemes show a better performance than the Arakawa and ex- 
plicit difference schemes. This clearly shows that designing a 
scheme with reduced truncation error or optimizing the scheme 
with respect to the dispersion relation could result in a better 
prediction capability than the designing a scheme with its dis- 
crete conservation properties. 
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Figure 24: Comparison of the second-order vorticity structure functions for 
Re = 1000 with a resolution of 1024 2 (Re c = 6.13) at time t = 5. 

In order to analyze the effects of the Reynolds number and 
required resolution to the performance of schemes, it is useful 
to define a cell Reynolds number in the following form 
2tt 

Re c =Re— (53) 

Therefore, in our previous experiments having a resolution of 
1024 2 at Re = 1000, the cell Reynolds number is Re c = 6.13. 



Table 7: Accuracy and efficiency of the finite difference approximations for the 
decaying turbulence problem at Re = 10 3 using a resolution of 1024 2 . The 
reference solution for the L2 norms is obtained by the pseudo-specral method 
which takes a CPU time of 621 .32 hrs. The speed-up ratio is defined as the ratio 
of CPU times. 



Scheme 


1Mb 


CPU (hrs.) 


Speed-up 


ED6 


1.96E-3 


128.39 


4.84 


ED4 


3.51E-2 


124.12 


5.00 


ED2 


1.07 


121.63 


5.11 


CD6 


1.63E-4 


149.52 


4.16 


CD4 


6.23E-3 


142.79 


4.35 


A4 


3.29E-2 


148.75 


4.17 


A2 


1.11 


127.29 


4.88 


DRP4 


1.80E-2 


126.59 


4.90 




Figure 25: Comparison of the second-order vorticity structure functions for 
Re = 1000 with a resolution of 512 2 (Re c = 12.27) at time t = 10. 



Next, we perform similar experiments for Re = 1000 with 
512 2 resolutions with an increasing cell Reynolds number of 
Re c = 12.27. The vorticity fields at time t = 10 obtained by 



nine different spatial schemes are plotted in Fig. 20 The cor- 
responding second-order vorticity structure functions are also 



illustrated in Fig. [25] We can easily see that the difference be- 
tween high-order and low-order schemes becomes amplified, 
and there is a substantial difference in the corresponding flow 
fields. If we decrease the resolution to the 25 6 2 resulting in a 
higher cell Reynolds number, Re c = 24.54, as shown in Fig. 21 
and Fig. [26] low-order accurate numerical schemes produce to- 
tally different flow fields. As the result of the series of calcu- 
lations performed in the second and third sets of experiments, 
we can see that the sixth-order accurate schemes always pre- 
dict the same accuracy results as the pseudospectral methods. 
In terms of fourth-order accurate schemes, we can see that the 
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compact difference scheme produces a slightly better accuracy 
than the dispersion-relation-preserving scheme due to its lower 
truncation error properties. Between the compact scheme and 
Arakawa scheme, the compact scheme exhibits better accuracy 
than the Arakawa scheme. The computational efficiencies in 
terms of speed-up ratios for these sets of experiments are very 
similar to the previous high resolution computations given by 
Table [7] Making a comparison between low-order and high- 
order schemes, we can see that a slight increase in the compu- 
tational time for a higher-order accurate scheme can result in 
significantly more accurate results than those obtained by the 
conventional second-order accurate schemes, and provide the 
same level of accuracy with the pseudospectral method in much 
less computational time. 




Figure 26: Comparison of the second-order vorticity structure functions for 
Re = 1000 with a resolution of 256 2 (Re c = 24.54) at time t = 6. 

We perform one more set of numerical experiments in order 
to demonstrate that the results of the various spatial schemes 
converge the same flow field under well resolved conditions 
using a resolution of 51 2 2 at a reduced Re = 100 for which 
Re c = 1.23. The vorticity contour plots obtained by these nine 
schemes are demonstrated in Fig. [22] which shows that there is 
negligible difference among them. Both low-order and high- 
order difference schemes predict the same flow field with the 
pseudospectral method. This can also be seen clearly from 
Fig. [27] showing the centerline vorticity distributions. Since the 
cell Reynolds number is smaller than 2, both second-order ac- 
curate and high-order accurate schemes produce similar results 
(i.e., well resolved direct numerical simulations). However, we 
can see that for larger cell Reynolds numbers, higher-order ac- 
curate schemes are required. 

6.4. Pile-up phenomenon for high cell Reynolds number 

Finally, we would like to test the performance of the finite 
difference formulations for extremely high cell Reynolds num- 




Figure 27: Comparison of the numerical schemes for Re - 100 with a resolu- 
tion of 512 2 (Re c = 1.23) at time t = 10. The centerline vorticity distributions 
at x = n are plotted. Inset: close-up of boxed area. 



bers. Two sets of numerical experiments are performed for a 
resolution of 256 2 at Re = 3000 and Re = 6000. The corre- 
sponding cell Reynolds numbers in these sets are Re c = 73.63 
and Re c = 147.26, respectively. A comparison of the fourth- 



order difference schemes is shown in Fig. 28 illustrating the 
evolution of the energy spectra for these Reynolds number. In- 
stead of presenting the angle averaged energy spectra, we prefer 



to show the full energy spectra, as defined by Eq. ( |50| , for each 
pair k = (k x , k y ) in wave space to make the comparison more 
visible. This figure clearly illustrates that a higher degree of 
pile-up phenomenon occurs for all the finite difference schemes 
except the Arakawa scheme. Some sort of aliasing errors pro- 
duce the growth of unphysical small scales. We demonsrate 
that these peaks of energy at small scales can be eliminated us- 
ing the conservative Arakawa scheme. Although there is no big 
difference between the Arakawa scheme and the explicit differ- 
ence scheme for smaller cell Reynolds numbers as shown in the 
previous sections, the discrete global conservation properties of 
the Arakawa schemes prevent the pile-up phenomena for higher 
cell Reynolds numbers. 

Another sets of experiments using a higher resolution of 
1024 2 are also performed to demonstrate the limits of the 
difference schemes by considering the pile-up phenomenon. 
Evolution of the energy spectra for both sixth-order compact 
and Arakawa schemes are shown in Fig. [29] Based on our 



numerical experiments, it is concluded that the pile-up phe- 
nomenon occurs with non fully conservative schemes for higher 
cell Reynolds numbers (Re c > 30). Therefore, the Arakawa 
scheme is a better candidate for underresolved simulations. It is 
mandatory to conserve energy to have a realistic simulations for 
higher cell Reynolds number. High-order centered schemes of- 
ten show spurious wiggles for convective dominated problems 
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Figure 28: Evolution of energy spectra in decaying turbulence for two sets of Reynolds numbers: Re = 3000 (Re c = 73.63) in the left column and Re = 6000 
(Re c = 147.26) in the right column. Comparison of the fourth-order finite difference schemes using a resolution of 256 2 : (a-b) explicit difference (ED4), (c-d) 
compact difference (CD4), (e-f) Arakawa (A4), and (g-h) dispersion-relation-pre^exving (DRP4) schemes. Energy spectra, defined by Eq. ( |50) , are shown for times 
t = 0, 2, 4, 6, 8, and 10. 
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Figure 29: Evolution of energy spectra in decaying turbulence using a resolution of 1024 2 obtained by (a) the sixth-order compact difference (CD6) scheme at 
Re = 1000 (Re c = 6.13), (b) the sixth-order compact difference (CD6) scheme at Re = 6000 (Re c = 36.78), (c) the forth-order Arakawa (A4) scheme at Re = 1000 
(Re c = 6.13), and (d) the forth-order Arakawa (A4) scheme at Re = 6000 (Re c = 36.78). Energy spectra, defined by Eq. {50) , are shown for times t = 0, 2, 4, 6, 8, 
and 10. 



when the numerical resolution is insufficient. The cell Reynolds 
number is restricted in some degree for these schemes. As the 
cell Reynolds number becomes larger, the nature of the finite 
difference equation is changed because of the cell size used. 
The spurious Nyquist signal would be generated almost instan- 
taneously that propagates into the whole domain resulting an 
instability in the solution. In order to prevent these numerical 
instabilities for a centered difference scheme without enough 
resolution, the filtering methodology can be also used for un- 
derresolved flows. 

7. Summary and conclusions 

A systematic comparison of a variety of high-order accu- 
rate finite difference schemes has been performed for two- 
dimensional decaying turbulence simulations by solving the 
vorticity-streamfunction formulation of two-dimensional in- 
compressible Navier-Stokes equations. The following schemes 
were considered: the 6th-, 4th-, and 2nd-order explicit dif- 
ference schemes, the 6th- and 4th-order compact difference 



schemes, the 4th- and 2nd-order Arakawa schemes, and the 4th- 
order dispersion-relation-preserving scheme. The objective of 
this study was to determine the accuracy and efficiency of these 
schemes for long-time integration of decaying turbulence simu- 
lations, which serves as an instance of a complex flow. We com- 
pared the schemes with the spectrally accurate Fourier-Galerkin 
pseudospectral method. A brief description of these algorithms 
is given in Section [3] In order to eliminate the errors coming 
from boundary conditions, we performed our experiments in a 
doubly periodic domain. This setting provides a suitable frame- 
work for testing the spatial discretization schemes. 

We first validated the methods by solving the Taylor-Green 
vortex problem which is one of the available exact solutions for 
two-dimensional unsteady incompressible flows. We showed 
that the theoretical orders of accuracies of the schemes were ob- 
tained in practice. Then we solved a more challenging double 
shear layer benchmark problem which has been used to test the 
grid-independence of many incompressible algorithms, since 
the presence of thinner and thinner shear layers evolving with 
time is not captured by low grid resolution representations. We 
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found that the low-order difference schemes require about twice 
the resolution in each direction to achieve the same results as 
the high-order difference schemes. We also demonstrated that 
numerical oscillations near the sharp vorticity gradients can be 
eliminated using the high-order difference formulations. 

We tested the efficiency and accuracy of several variants of 
third- and fourth-order Runge-Kutta time stepping algorithms 
for decaying two-dimensional turbulence. The numerical sta- 
bilities of the spatio-temporal schemes were addressed by per- 
forming a time refinement study. We first showed that the 
fourth-order Runge-Kutta schemes have greater accuracy for a 
given time step and have larger allowable stability regions than 
the third-order Runge-Kutta schemes. We found that the spa- 
tially second-order accurate schemes can run using almost triple 
the effective time step compared with the sixth-order schemes, 
and can also run using almost double the effective time step 
compared with the fourth-order schemes, thereby making them 
more competitive with the high-order finite difference schemes. 
On the other hand, we showed that the second-order schemes 
usually require double the resolution in each direction to be able 
to obtain a similar accuracy. Therefore, we found that higher 
order schemes become more effective in terms of the tradeoff 
between the accuracy and efficiency. 

We studied the Reynolds number dependence of freely evolv- 
ing decaying turbulence. We showed that the predicted energy 
spectrum asymptotically converged to the theoretical k~ 3 scal- 
ing as the Reynolds number increased, which is predicted by 
the KBL theory for forward cascading two-dimensional tur- 
bulence in the inviscid limit. We also demonstrated that the 
analysis was also relatively free from temporal discretization 
errors. We then tested the spatial finite difference schemes for 
various Reynolds numbers. We defined a cell Reynolds num- 
ber Re c = Re2n/N x . We demonstrated that for well resolved 
simulations, accurate results are obtained for all the schemes if 
Re c < 2. If we increase the Reynolds number, however, the 
difference between the high-order accurate and second-order 
accurate schemes increases dramatically. Compact difference 
schemes provide more accurate results if Re c < 30. Spurious 
Nyquist signals have been observed for larger cell Reynolds 
numbers resulting in a pile-up phenomenon. In this case, the 
fully conservative Arakawa schemes gives more accurate re- 
sults. 

Our results demonstrate the importance of high-order ac- 
curate representations in convection-dominated complex flow 
problems, which, in general, require long time integrations. 
Although the order of accuracy is less important for well- 
resolved direct numerical simulations, it becomes significant 
for Re c > 2. We show that, for fully resolved simulations, 
high-order accurate schemes do not require spectral-like sim- 
ulation times. In fact, spectral-like accuracy is obtained with 
a speed-up factor of 5 over the pseudospectral method using 
the 6th-order compact scheme. We also demonstrate that the 
fourth-order Arakawa scheme is a better choice for higher cell 
Reynolds number computations. 
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